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Abstract 

Some statistical properties of finite-time stability exponents in the standard map can be es- 
timated analytically. The mean exponent averaged over the entire phase space behaves quite 
differently from all the other cumulants. Whereas the mean carries information about the strength 
of the interaction, and only indirect information about dynamical correlations, the higher cumu- 
lants carry information about dynamical correlations and essentially no information about the 
interaction strength. In particular, the variance and higher cumulants of the exponent are very 
sensitive to dynamical correlations and easily detect the presence of very small islands of regular 
motion via their anomalous time-scalings. The average of the stability matrix' inverse trace is even 
more sensitive to the presence of small islands and has a seemingly fractal behavior in the standard 
map parameter. The usual accelerator modes and the small islands created through double saddle 
node bifurcations, which come halfway between the positions in interaction strength of the usual 
accelerator modes, are clearly visible in the variance, whose time scaling is capable of detecting 
the presence of islands as small as 0.01% of the phase space. We study these quantities with a 
local approximation to the trace of the stability matrix which significantly simplifies the numerical 
calculations as well as allows for generalization of these methods to higher dimensions. We also 
discuss the nature of this local approximation in some detail. 
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I. INTRODUCTION 



A central and widely used measure of instability in dynamical systems is the Lyapunov 
exponent . If a positive Lyapunov exponent is found for a continuous family of trajecto- 
ries, then it is taken as a measure of their dynamical chaos. It is natural that the Lyapunov 
exponent, which is defined as a time-approaching-infinity limiting value, when calculated 
numerically, would be approximated by considering finite time segments of the infinite his- 
tory of a trajectory. This has led to the study of the so-called finite time or local Lyapunov 
exponents (FTLE) [2I in a variety of contexts. For example, FTLE in one-dimensional maps 



0, Q, y, 0], Hamiltonian systems [?, y, [q], advection and turbulent flows 1^, 
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just a part of a rather large literature. Additionally in most experimental situations one 
encounters finite time quantities. Such FTLE have been the subject of numerous studies 
and have been experimentally measured in many contexts such as for instance in the study 
of passively advected particles p^]. There is also a very intimate relationship between the 
theory, methods, and measures found in the subject of Anderson localization [13| and FTLE. 
Although, typically dynamical correlations are assumed vanishing or if correlations are ac- 



counted for 



14| . then they are not the kind which is found in Hamiltonian flows where all 



the complexity of boundaries between chaotic and regular motion may exist. In other words 
these correlations are then calculated assuming that a typical trajectory covers phase space 
uniformly, whereas the presence of islands, however small, ensures that this is not true. 

A very closely related, but less directly studied, quantity is the stabihty exponent which is 
defined via the eigenvalues of the stability matrix. A detailed study of this measure vis-a-vis 
the Lyapunov exponent was carried out in [15] where it was conjectured that the Lyapunov 
and stability exponents were equal under a broad range of dynamical assumptions. It turns 
out that some of the studies that purported to be about the FTLE have actually been studies 
of the stability exponents instead. In semiclassical mechanics, such as in the Gutzwiller trace 
formula [igI], one of the natural classical quantities that appears is, in fact, the stability 
exponent and not the Lyapunov exponent. In higher dimensional systems (larger than two- 
degrees-of- freedom) , it is the sum of the positive stability exponents that appears and in the 
infinite-time limit this tends to the Kolmogorov- Sinai (K-S) entropy. Since it has long been 
appreciated that the small wavelength limit and the large time limit in quantum mechanics 
are not interchangeable, we may expect that for fixed wavelength, it is the finite-time stability 
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exponents, to be defined below more precisely, which may be relevant; see for example [iTJ. 
In fact, many of the interesting ways in which the fluctuations in the stability exponents 
might potentially show up in quantum mechanics and semiclassical theories remains to be 
investigated. This forms our primary motivation for studying the stability exponents, their 
distribution and their scaling properties. Nevertheless, in this paper we restrict our attention 
to purely classical and low-dimensional systems (two-dimensional area preserving maps) to 
establish basic properties and approximations, as well as to demonstrate their usefulness 
in detecting small stable dynamical structures. Naturally there are many parallels in the 
large literature on FTLE both of a quantitative and qualitative character given that for a 
broad class of dynamical systems the infinite-time stability exponents are identical to the 
Lyapunov exponent. Our results presented here generally for stability exponents are thus of 
relevance to the literature on the FTLE as well. 

In the following section we first introduce the model in use throughout (the standard map) 
and an approximation for the stability exponent which provides a practical, calculational 
simplification that can be generalized to higher dimensions. Its ergodic value has corrections 
to the well-known formula of Chirikov 18|]. In section UTTl the local properties of the stability 
exponents are studied. We image the finite-time stability exponents in phase space and 
show their close relationship to the stable and unstable manifolds of the map, and show 
their connection to locating small non-hyperbolic regions. In section llVt the fluctuations 
of the stability exponents are investigated, and a formula for the ergodic average of the 
variance is derived. The scaling of the variance is shown to be a sensitive sensor for small 
islands. It is also shown that quantities that incorporate higher cumulants and that are of 
natural semiclassical relevance are even better sensors of such phase-space structures. In 
the appendices, formulas are collected for an expansion of the exact trace of the stability 
matrix as well as results for the time scaling of the stability exponent higher cumulants. 



II. BACKGROUND, AN APPROXIMATION, AND ERGODIC AVERAGE 

In this section we introduce the kicked rotor, provide definitions of the quantities of inter- 
est, explore an approximation of great practical utility, and use it to generate an improved 
analytic estimate of the Lyapunov exponent for the standard map. 



3 



A. The Hamiltonian 

Consider the 2-D standard map well-suited vehicle for our investigations. 

The methods and conclusions we reach will have wider applicability as this map represents 
the generic, local dynamics of any integrable system with shear in action-angle coordinates 
and shows how that structure gets altered by perturbation. A general kicked rotor is a 
mechanical-type particle constrained to move on a ring that is kicked instantaneously every 
multiple of a unit time, t = nr. Supposing the radius of the ring to be l/27r and r = 1, the 
Hamiltonian takes the form 

H{q,p) = ^ + V{q) 5{t-n). (1) 

n=— oo 

where V{q) is a function periodic on the interval q G [0, 1). From H{q,p), mapping equations 
can be given: 

Pi+i = Pi- V'iqi) mod 1 

qi+i = qi+pi+i modi. (2) 

where we have made the choice to apply the potential kick before the free motion. The 
notation V indicates the derivative of V with respect to q. We have also restricted the 
phase space to the unit torus {q,p G [0, 1)) for convenience; this is of no consequence for the 
results presented in this paper. The simplest periodic function on a ring is just the lowest 
harmonic and leads to the standard map. 

V{q) = -^cosi2nq) (3) 



Many results have long been known for the standard mapjlSj. At = the map is 
integrable and is essentially a stroboscopic map of a freely rotating particle. There are both 
rational and irrational tori, depending on whether the frequency of rotation is commensurate 
or not with the frequency of the strobe. As the kicking strength is increased from zero with 
exactly the same frequency of stroboscopic observation, the incommensurate tori (irrational) 
survive small perturbations in accordance with the KAM theorem, while the commensurate 
(rational) ones break up into a pair of stable and unstable orbits in accordance with the 
Poincare -Birkhoff theorem. The phase space becomes mixed with stable and chaotic orbits 
for increasing K. At around K ^ 1 the last rotational irrational KAM tori breaks and 



this leads to global diffusion. Up to = 4 a stable fixed point persists. Beyond K ^ 5 
the standard map is considered to be largely chaotic, although it is also not proven to be 
completely chaotic for any value of K as far as we know. At an infinity of values of K 
stable fixed points are known to appear in the p = line for the map on the torus that are 
accelerator modes fa the map on the p,hnder fl. These typically occupy tegrons m the 
phase space whose areas scale as [18]. Otherwise, the Lyapunov exponent of the map, 

to a very good approximation (Chirikov's result), increases with K as ln{K/2). 



B. The finite time stability exponents 

For the map of Eq. ([2]), the stability matrix at time t is given by 

t—i / 1 V^"f I 

to to \ l l-V"{qi)J 

where the M{qi) are local one-step stability matrices. The finite time Lyapunov exponents 
are defined via the singular values of Mf, which are the eigenvalues of Mj Mt and turn out to 
be non-negative real numbers. The eigenvalues of stability matrices Mt on the other hand can 
in general be complex. To be technically correct we refer to the positive finite-time "stability 
exponents" (FTSE), \t{qo-,Po)-, where the time t is the integer number of applications of the 
mapping equations and {qo.Vo) is the initial condition. Consider only initial conditions 
{qo,Po) for which the orbits are unstable, in which case for 2D area preserving maps such as 
we are presently studying the eigenvalues have to be real. The Xt{qo,Po) can be expressed 
in terms of the trace of the stability matrix as follows 

iln|Tr[Mi(go,Po)]| (5) 

where Tr(...) denotes the trace operation. The approximation applies for t large enough. 
In that case, the approximation error is exponentially small in the time except for stable 
trajectories. The generalization to higher dimensions is through the exponential growth rate 
of the magnitudes of the eigenvalues of the stability matrix (the Jacobian) at time t. 

The natural quantity, which occurs in semiclassical expressions, is |Det(M( — /)| whose 
large time behavior is determined by the sum of the positive stability exponents. For 2D 
area-preserving maps one has |Det(/ — Mt)\ = |Tr(Mj) — 2|, and the maximum stability 
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K{qo,Po) = 



Tr[Mi(go,Po)] + VTr[Mi(go,PoF-4 



exponent coincides with the K-S entropy in the t — oo hmit. Thus for the kicked rotor, 
studying |Det(/ — Mt)\, |Tr(Mt)|, or the K-S entropy is equivalent and our interest is in 
these quantities' (or their logarithms) finite-time fluctuations over phase space. Therefore 
our study is naturally close to the well-explored fleld of FTLE. There are strong reasons to 
believe that in the inflnite time limit the deflnitions of FTLE and FTSE are equivalent in 
the chaotic cases that we are concerned with, although non-generic examples are known for 



which this is not true 



15| . Since we are dealing with finite time measures we wish to be 



technically correct and refer to FTSE, although these typically deviate from the FTLE by 
terms of the order of l/t. Recall that the FTLE can also be defined via the way initial error 
vectors yo grow in length: 

/^H?o,Po;2/o) = Tin (6) 

V w yo \\ J 

This definition of FTLE has the added complication that they will depend on the initial 
error vector yo as well. However, in the infinite time limit this corresponds to the exponent 
coming from the singular values of the limit of Mt-,oo- Thus, fi^o is the Lyapunov exponent 
and it turns out not to have initial condition dependence or fiuctuations. There are reliable 
methods of evaluating each of these quantities [l^ (whether FTLE or FTSE) and they 
give essentially the same information, especially regarding the fiuctuations of the exponents 
that we are primarily concerned with here. Higher dimensional generalizations, which are 
currently under study [2^, involve the fiuctuations of |Det(M — /)| and therefore we prefer 
to concentrate on FTSE in this paper on 2D maps. 



C. Approximating the trace 

Before turning to a more detailed analysis of Xt{qo,po), its average, variance and higher or- 
der moments and cumulants, we investigate approximations to the trace. These approxima- 
tions make the trace calculations both efficient and simple, and extend to higher-dimensional 
systems, where they provide significant advantages. For any stability matrix of the form of 
Eq. (HI), it turns out that Tr[M((go,Po)] can be expressed in an exact, closed form of products 
of the individual traces Tr[M(gj)] = 2 — V"{qi)] the full expressions to all orders are found 
in Appendix A. We display here the zeroth order approximation, indicated by j = 0, and 
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its first order correction, indicated by j = + 1: 

t-i 

Tr[Mi(go,Po)],=o = H - ^"(^0] (7) 

t—l / cyclic 

T.|M,(*.p.)l,.„,. = n 12 - ^"M] X - E |2-t^»fa)||2-V»fe.)| 

These expressions represent well the large- i^' behavior of the exact trace TT{Mt{qo,Po)), even 
though, a priori, one might suspect that estimating this trace by Eq. is rather inadequate. 
After all, consider the fact that the leading correction, shown in Eq. ([8]), contains t terms 
of inverse products of pseudo-random factors which have the possibility of nearly vanishing. 
Crudely speaking, as multiple trajectories are considered or as t increases for an individual 
trajectory, it becomes more and more likely that a nearly vanishing factor will be encountered 
thus creating a correction term that dominates the leading term. Worse, all the successive 
corrections involve combinations of the same products of inverse factors (Appendix A). It is 
possible that no particular correction term dominates either. 

Fortunately and oddly enough, Eq. ([7]) (the j = term) is an excellent approximation, and 
in fact, is more than adequate for analytically estimating the leading asymptotic behaviors 
of the finite time stability exponents. It is even good enough to capture the power-laws for 
both local averages and variances of FTSE, which we investigate ahead. Figure ([I]) compares 
for a single trajectory the natural logarithm of the exact trace, the product form estimation 
of Eq. ([7]), and the estimation of Eq. ([8]) where the mean exponential increase has been 
subtracted from all three in the first panel. In these panels several general features are 
readily seen. As a function of time, the leading approximation follows the exact trace time 
step by time step except for sudden vertical displacements at time steps where the factor 
|2 — K cos{2'Kq)\ nearly vanishes. In this example, the corrected approximation follows the 
exact trace longer in time before departing; although, this does not have to be the case. A 
second point is that the exponential growth dwarfs the errors of the approximation. Even 
though the approximation of Eq. ([7j) is off by a factor of roughly ~ e^^^ before t = 300, this 
is miniscule compared with the average exponential growth of the trace (in this example it 
is exp [t In ii'/2] = 4*) and so it alters the stability exponent very little. From the form of 
the leading correction, a Levy-flight nature might be expected for the corrections. This can 
be seen in the natural logarithm of the ratio of Eqs. (17f8|) respectively where only a couple 
of jumps dominate the differences between the zeroth and first approximation. Although 

7 





FIG. 1: Comparison of the two approximations, Eqs. (j7l8p . for the trace Mt{qo,po) with the exact 
value. In panel (a), the three curves are for K=8 with the initial condition {qo,Po) = (0.3,0.6), 
A = \n\Mt{qo,po)\ — tln(K/2), the exact trace is the solid line, the j = approximation is 
the dotted line (lowest), and the j = + 1 approximation is the dashed line. In panel (b), r = 
In \Mt{qo,Po)\/\MtiqoiPo)\ ■ For the solid line, Mtiqo,Po) is evaluated to the j = approximation, 
and for the dashed line, Mt{qo,po) is evaluated to the j = + 1 approximation. Similarly to panel 
(b), in panel (d), the ratio of the two approximations is given. Panel (c) is similar to panel 
(b), except that K = 9.28 with the initial condition {qo,po) = (0.19,0.38) (which is a stable 
trajectory). The solid line is for the j = approximation whereas the lower curve is for the 
j = + 1 approximation. 

not shown, as K increases, a corresponding curve seems to be rather "quiet" punctuated by 
relatively fewer jumps. 

These approximations naturally are worst for segments of trajectories that are stable or 
nearly stable. For them 2 — V" is repeatedly close to zero. The large value of K thus does 
not order corrections. An example is given of a nearly stable trajectory in a region near a 
double saddle point bifurcation. The value of K is larger than used previously to illustrate 
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the behavior of a highly unstable orbit. Nevertheless, the approximation is far worse. On 
the other hand however, it gives a stability (or Lyapunov) exponent of about 0.067, a small 
number compared with fi^o ~ 1.53. As long as the percentage of stable and nearly-stable 
orbits is very small, the errors in using Eq. ([7]) for all the orbits regardless of their nature 
are small after averaging over initial conditions. For the usual Lyapunov exponent there 
exists a "local" expression which is in some sense the equivalent of the j = expression for 
the trace we have just examined ^j. 



D. The average finite-time stability exponents 

Much of what is known applies readily if a product form can be found for Tr[Mj(go,Po)]- 
This is precisely what the kicked rotor has for sufficiently strong potentials and times not 
too long. There the j = term would be sufficient, and from Eqs. fl5|7|) . the {At(go)Po)} are 
approximately given by 

A?(go,Po) = ^£ln|2-r"(g.)| (9) 

i=0 

The superscript indicates that it is the zeroth order approximation. Assuming that the 
system is fully chaotic, the expectation value is given by integrating (averaging) over all 
initial conditions and 

(A?)= [\q ['dpX',{q,p) = lYl rdq\n\2-V"iq)\= f dq \n\2 - V" {q)\ (10) 

JQ Jo "^0 



For the standard map, the result is (A*^) = In^; the integral is evaluated exactly. Not 



surprismg. 



: ly, this coincides exactly with the Lyapunov exponent approximation made by 



Chirikov 18| even though this result is for the stability exponent. The phase-space mean 
of the finite-time stability exponent is independent of integration time and time is dropped 
from its notation. This indicates that the mean is to within approximations made so far 
equal to fi^o- In contrast, the entire phase space average of the FTLE is not equal to its 
infinite time ergodic average, but has a transient that decays as a power law in time and in 
the absence of significant islands of regular motion goes as 1/t. 

We plot the mean stability exponent in Fig. ([2]). It is interesting that there is a simple 
systematic correction to the well-known formula of Chirikov above that is of the order of 
l/K"^, and that other than the value of K (excluding sample size noise), there is essentially 
no indication of any other kind of dynamical information. 
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FIG. 2: The exact mean stability exponent calculated as a function of kicking strength with the 
zeroth order approximation, ln(K/2), subtracted. Other than small oscillations, the difference is 
seen to be almost exactly equal to the first correction, (X^ — 4)~^, given in Eq. (jl6l) . This is plotted 
as the dashed curve . 

The origin of the corrections lie in the terms with j 7^ in the approximations to the exact 
trace, however its rigorous derivation seems fraught with difficulties. A somewhat heuristic 
derivation of the l/i^'^-correction to Chirikov's expression can be derived from Eq. ([8]). The 
product term gives an independent contribution and is exactly lnK/2 as already noted. 
Subtracting that term, one is left with the ergodic average of at least the j = 1 term to 
determine the correction; 

^ / cyclic ^ 

SX{K) - ^ 1 - XI cos 27iqi) {2 - K cos 2nqi+i) 

We are not aware of a way to evaluate this integral without making a few simplifications. 
Consider first that the correction is essentially being expanded about = 00 or in other 
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words the expansion is about e = l/i^ = 0. In this hmit at some fixed value t, the probabihty 
is vanishing that more than one factor of |2 — i^' cos 27iqi\ for a given trajectory is close enough 
to zero that it can meaningfully contribute. Let go denote the q value for the minimum of 
the set of |2 — K cos 27rgj | . Although that involves two terms in the series, only one will 
dominate. However, since there are t-terms in the cyclic sum, there are t ways this can 
happen; together this is consistent with what happens in the subject of extreme statistics. 
The basic expression to be evaluated is 

1 



5X{K) 



In 



1 - 



(12) 



(2 - Kcos27Tqo) (2 - ii'cos27rgi^ 
Due to the momentum not entering the expression, {qo, qi} can be integrated as two inde- 
pendent variables each uniformly distributed on the interval [0,1). The next step follows 
from the discontinuity in the expression for the integral 



dgln |a ± 6 cos 27rg| 



In 



\b\ 



In 



(13) 



\b\ > \a\ 
2 l«l > 1^1 

It is straightforward to verify that to the extent that the upper form of the integral applies, 
the expression for 6X{K) vanishes. Therefore it turns out that only values of go for which 



< 2 — K cos 27rgo < 



(14) 



K -2 - - K + 2 

contribute a correction to the average stability exponent and the limits of integration on go 
are restricted to that domain (actually there are two such domains). Switching variables to 
y = 2 — K cos27iqo, one has (accounting for both domains) 



6X{K) 



71 



K+2 



^ ^K^ - (2 - yf 



In 



l-2y+ v/(l-2y)2-irV 



In 



K\y\ 



(15) 



If the leading square root factor is expanded in a series of powers of y, the integrals can be 
analytically evaluated. This gives, 

1 19 4 24 



5\{K) 



+ 



+ 



+ 



+ 



(16) 



K2-4 6(7^2-4)3 (i^2_4)4 (i^2_4)5 

As seen in Fig. the first term accounts completely for the overall correction to the 
Chirikov result to the accuracy level of the numerical calculations. The j = 1 trace correction 
appears to have given properly the stability (Lyapunov) exponent to 0{K^^). This is 
interesting because it seems that it is not necessary to consider the j = 2 and higher 
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corrections to the trace whereas one might have naively expected the j — 2 trace correction 
to give an 0{K~^) contribution. Chirikov's accurate estimate for the Lyapunov exponent 
also has the same corrections derived above. This may be expected as at infinite times 

the Lyapunov and stability exponents indeed converge for the standard map. However, 
numerical calculations of the FTLE, even averaged over the entire phase space, mask this 
with transients that are of the order of 1/ {tK) . 

III. LOCAL PROPERTIES OF FINITE-TIME STABILITY EXPONENTS 

The FTSE Xt{qo:Po) is also a local stability exponent as emphasized by its dependence 
on initial conditions. The mean over a local "area" A that is much smaller than the entire 

phase-space (Q) includes transients that arc power-laws in time. Its interesting variations 
over phase space can be displayed by assuming the following particular form: 

{MQo,Po))a = ^ + HK/2) + + 0{K-% (17) 

The area A is ideally a very small local region of the available phase space centered at (go, Po) 
and the final two terms are responsible for removing yUoo, thus leaving just the fiuctuating 
part. Here K » 5 (in practice K > Q) and the time t is larger than the logarithmically 
short mixing time scale, = \n{Vt/ A) / Hoo ~ ln(l/^)/ln(i^/2). This time scale arises as 
the "log-time" during which the initial region of volume A through propagation essentially 
elongates to the extent that it enters each distinct local region of area A once on average. 
The averaging in the LHS of the above expression is over the uniform Lebesgue measure, 
which is the invariant measure of the map. Note that when the entire phase-space, i.e. the 
unit square for the kicked rotor, is used for averaging the FTSE, the transients in time seem 
to be exponentially small and essentially Cq — 0. 

Expressing the finite time stability exponents in this way allows us also to draw attention 
to anomalous features that may be present for various values of especially those with 
small, but significant, stable islands in phase space that survive for large such as the 
accelerator modes. The quantity Cjs, fluctuates both in time and phase space point around 
which the averaging is done, it can also have a strong dependency on K. The index 7 serves 
to remove the principal non-fluctuating time scaling part of the FTSE and determines the 
rate at which it approaches the ergodic average. For example, in the case when there are 
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FIG. 3: A density plot of the average (A-i-io(q'o,Po))^ (top, bottom respectively) as a function of 
(go,Po)- High to low density is indicated as a color transition from purple to red. The area A over 
which the averages are calculated are squares of area 0.0052. rp^g 

case = 10 (center column) 

has no apparent stable islands, whereas at = 9.26 (left column) there are two small islands each 
of period four where the density is enhanced and purple (see the next figure for a magnification of 
the dynamical structure underlying the purple zone). Note the prominent linear structures for a 
large value of K (=40) in the right column. 

no apparent small islands 7 = 1 (see text ahead), while in the presence of stable regions 
< 7 < 1. The area A serves as a region of coarse-graining and it must be emphasized that 
unless we choose A to be the entire phase space (fi), the regions will exclude regular islands 
and trajectories. We present results for an exponent similar to 7 when studying the variance 



below. Previous studies of the FT 
infinite time Lyapunov exponents 



E have a. 
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so presented such forms of convergence to the 



22| . but typically in dissipative systems. 



To get an idea of the complexity of the phase space distribution of the FTSE as we 
move (go, Po) around, always averaging over a small area A of phase space centered at this 
point, we display density plots in Fig. (I3l). The figures show forward and backward in time 
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evolution for three cases of namely 9.26, 10, and 40 respectively. The intricate structures 
that arise highlight areas of reduced hyperbolicity and hence smaller FTSE. For the case of 
forward evolution in time, trajectories which are trapped in more stable areas or propagate 
into those areas have lowest stability and appear purple. Those trajectories propagating 
into more stable zones must do so in a direction aligned along the greatest exponential 
compression. This direction is determined by the stable manifolds of the short periodic 
orbits. Thus, one sees a structure that mirrors the foliation of phase space by these stable 
manifolds in the top row of the figure. Similarly, for reverse evolution in time, one sees a 
structure that mirrors the corresponding unstable manifolds in the lower row. In the first 
column, the phase space has small islands, while the center column, = 10, does not have 
such prominent islands. The absence of a detectable island of regular motion results in a less 
contrasted variation of the FTSE. At = 40 (right column), which is very strongly chaotic, 
there are linear regions of low FTSE that are clearly visible. However, any other regions 
have more or less merged into the background. These linear regions coincide with the K 
independent borders of the non-hyperbolic regions for a one-step forward (the vertical lines 
at g = 1/4 and q = 3/4) and one-step backward iteration (the lines q — p — 1/4 = Omodl 
and g — p — 3/4 = mod 1). The fact that these lines are not dense in the phase space and 
survive for long times at large K underlies the coexistence of strongly chaotic dynamics with 
the fact that the map may not possess a single K value for which it becomes mathematically 
fully hyperbolic. Similar plots for the FTLE have been plotted earlier in the contexts of 
area-preserving maps for smaller values of K and passive scalar advection 11|. 

It is possible to use this plot to help locate small islands. Investigating the purple zones 
seen in the left column of Fig. ([3]) more closely reveals the interesting structure found in 
Fig. (jl]). It arises for K values greater than approximately 9.2 at which point there is a 
double saddle node bifurcation. As K increases from the bifurcation, the structure grows 
to a maximum size near K = 9.26. The full structure includes not only the small stable 
islands, but also the short unstable periodic orbits along with their homoclinic tangles, i.e. 
crossings of stable and unstable manifolds. The turnstiles formed by these crossings limit 
the rate trajectories can move in and out of the region surrounding the islands of stable 
motion (24^. Therefore in this structure, it is not just the immediate interface between 
the chaotic and regular motion which is generating trapping, although that, of course, is 
happening additionally just beyond the outermost tori pictured. In fact, chaotic trajectories 
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0.18 0.20 0.22 0.24 0.26 

q 

FIG. 4: Expanded plot of two islands of period- four regular motion and the stable and unstable 
manifolds associated with two unstable orbits born at the same double saddle node bifurcation 
(K=9. 20667769750) as the islands themselves. Only one of this structure's four images is plotted. 
Transport in and out of the inner region is severely limited by the small turnstile fluxes visible. 

may remain trapped inside the tangle for hundreds or more time steps in the roughly 2/3 
portion outside the islands themselves. For this reason, the entire structure was visible in 
Fig. ([3]), not just the islands, and multiple trapping time scales would be expected. One of 
the more unusual features is a highly asymmetric instability (factor of 1000 or so) in two 
of the stable and unstable manifolds pictured. The interior middle line is actually a stable 
and unstable manifold whose turnstile is so tiny as to appear as a single line dividing the 
interior, thus making transport between the upper and lower region nearly non-existent. 



15 



IV. FLUCTUATIONS OF THE FINITE-TIME STABILITY EXPONENT 



The statistical properties of {At(go,Po)} can be studied via the construction of a proba- 
bihty density pxXx'i i-e. the probabihty of finding in the interval dx centered at x. For 
a large class of chaotic systems, a great deal is known about general properties of px^{x;t), 
although strictly speaking just for the FTLE. The approach to a unique Lyapunov exponent 
is connected to p\t{x; t) approaching a 5-function limit as t ^ oo. Through the study of the 
density's cumulants for fully hyperbolic systems, one sees that there is a slow approach to 
an ever narrower, roughly Gaussian density as t — > cxo [2I. A detailed study with analytic 
results has also been performed for "white- noise" random media 25|]. Interestingly enough, 
it is also known that this limiting process is too slow to understand the behavior of the 
generating function for the cumulants z/„(t) 

In(expmtAi) = f^z/„(t)^^ ^ ln(|Tr(Mi)r) (18) 

n=l 

which is also a quantity of direct interest in semiclassical studies or the theory of Anderson 
localization. This relation should be thought of as an expansion in the real variable m. 

If one were to calculate the full probability density px^ {x; t) for the kicked rotor, one would 
see that it qualitatively takes on one of two forms. If the kicking strength is set to a value at 
which there is no discernible evidence for any stable orbits (no detectable non-exponentially 
decaying dynamical correlations), the density has a well defined peak near the Lyapunov 
exponent, and is not too far off a Gaussian form. If the kicking strength is set to a value 
at which both stable and strongly chaotic trajectories co-exist in the phase space, then it 
takes on a bi-modal form with a peak near the mean Lyapunov or stability exponent and 
another representing algebraic behavior near zero Qj. The density is continuous between 
the two peaks. Nevertheless, as t — > 00, a unique Lyapunov exponent emerges albeit much 
more slowly. This has a profound effect on the time scaling properties of the density or more 
specifically its cumulants (moments) with the exception of its mean. The algebraic peak 
and the 'neck' connecting the two peaks disappears in the t ^ 00 limit very differently than 
the rate at which the width of the main peak narrows. Time scaling therefore gives a great 
deal of information about non-hyperbolicities. 
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A. Time scaling of the variance 



We proceed next to a more detailed analysis of the variance (second cumulant) by first 
deriving the contributions that neglect dynamical correlations, which leads to an explicit 
formula for the standard map; the first cumulant (the mean) is already treated above. It is 
shown that while the general time scaling of the variance is this is significantly modified 
by the presence of even very small stable islands, thus providing a practical way of detecting 
a 'barely' mixed phase space. We anticipate that this would be especially useful in its 
ligher-dimensional generalizations [2o|. A first study appears to support our expectations 



26|. 



The variance of the FTSE over Q can be expressed simply using the approximation of 
Eq. ([7j) for the exact trace; the exact one, z/2(t), differs from its approximation above, 1^); 
by very small amounts which can be ignored. One finds: 

^2it) = dqodpo (^7 iiln |2 - V"{qi)\^ - {v[f (19) 

The variance may be be separated into a term that is explicitly independent of dynamical 
correlations and a series of terms that contain the contributions of dynamical correlations, 
which respectively are: 



i,l(t) = i / / dpd, (ln= |2 - V"{q)\ - [vlf) 
'0 Jo 

t-l „1 „1 



/ dpodgo(ln|2-V"(go)|ln|2-V"(g.)|-(z/°)2) (20) 
^ ^=1 "'0 Jo 

The contribution of the first term above is 

-X2^- / dg(ln|Kcos(27rg)|-ln(ir/2))2 = -^ = -0.822.... (21) 

Here we have introduced X2 as a term that does not carry correlation information. 

Depending on the system, strong behavior connected to dynamical correlations but not 
reflected in the mean, is possible in the variance, which the summation above makes clear. 
Note that the dynamical correlation function is bounded for a broad class of potentials. 
Furthermore, as long as the correlations decay faster than r~^, they cannot alter the basic 
time-dependence of the variance. If they decay more slowly, they could well become the 
leading contribution in time. As seen in Fig. ([5]) the variance is strongly dependent on the 
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FIG. 5: The j = order calculation of the variance as a function of kicking strength. It is 
multiplied by t and averaged over t = 150. It illustrates that the variance is very sensitive to 
dynamical correlations. The dashed line at 7r^/12 is the value of the integral from Eq. (j21|) . 

value of the kicking strength K and shows marked peaks at various values. These include 
conventional accelerator modes as well as other stable islands, such as those illustrated in 
Fig.(jl]) that appear near K = 9.3. 

B. The variance as a sensor for small elliptic islands 

As just noted, the variance scales with time differently from the 1/t law above as a result 
of correlations decaying as slow or more slowly than r^^. For the kicked rotor, consider 
ensembles of trajectories that originate from small regions A of the phase space rather than 
the whole space; the reduced averaging region is motivated by the interest in a process 

18 




that more practically generalizes to systems with greater numbers of degrees of freedom. 
The average behavior (A^(go,Po))y4 may obscure the scaling by the presence of the infinite- 
time asymptotic value, whereas the fluctuations die out asymptotically and are more stable 
for evaluating power laws. Although, note that "wiggles" in the FTLE were attributed to 
stability islands in an earlier work on kicked tops [27]. The variance can be written 

Var[A°(go,Po)]^ = /^^/t". (22) 

where (^) = Var[A°(go,Po)]Q- After accounting properly for the dependence, one 
sees that fluctuates as a function of time similarly to a diffusive variable. Generally 
speaking, the greater the time range over which a is calculated and the greater the number 
of trajectories run, the greater precision with which a is determined. In this regard, the 
simplicity of the j = approximation is quite helpful because it makes it possible to increase 
both in a practical sense. One suspects that the smaller and less influential the island of 
regular motion, the closer a approaches unity and thus a-accuracy is directly linked to how 



29| have studied 



small an island one can detect in this way. Earlier works such as |28l . 
scaling laws that arise in the fluctuations of the FTLE and have already emphasized the 
role of sticky regions and islands. Ahead, we show islands of measure roughly 0.01% of Q 
that were found with the FTSE time scalings relatively easily. 

In Fig. ([6]) we show the scaling of the FTSE variance for a few values of the kicking 
parameter compared with a = 1. The regular island structures of Fig.Q at i^' = 9.26 
have a very signiflcant effect on the time scaling. One also sees the appearance of more 
than one time scale, i.e. slope, depending on the time regime. The turnstiles previously 
pictured would at flrst lead to only small numbers of trajectories entering into the interior 
region of the homoclinic tangle, and one would expect a to begin close to unity. Later, a 
certain proportion would get trapped inside for up to hundreds or thousands of time steps. 
During this period, the greatest proportion would resemble nearly stable orbits and a would 
deviate most from unity. Further in time, only the trapping right at the regular island-chaos 
interface would contribute to deviations from unity and the slope increases somewhat again. 
It is known that different trapping mechanisms can co-exist in phase space and lead to a 
"multifractal" process jsol . This particular example is instructive because the turnstile time 
scale is fairly well separated from the others. The smaller the stable islands, the more the 
power laws recover the uncorrelated \jt form. KX, K = (8.47,9.35) there are only small 
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FIG. 6: The scaling of the variance for an ensemble of 40000 initial conditions uniformly distributed 
in the square A = {(0.10, 0.10), (0.15,0.15)}. 

detectable deviations from unity, and in both cases it is possible to locate very small islands; 
see Fig. ([71). At K = 9.5 we could not locate islands in this way and the exponent a is as 
close to unity as the precision of the calculation. 

It is worthwhile remarking that while accelerator modes give rise to anomalous diffusion 



(in momentum) in the standard map on the cylinder j3l|, the different scaling regimes 
discussed above hold whether the stable regions are accelerator modes or not. However, 
the suppression of the rate of the variance's approach to zero and the average Lyapunov 
or stability exponent to the ergodic average in the presence of small islands has similar 
origins. Namely, they are due to trajectories spending intermittently long intervals of time 
near sticky islands or regions of almost marginal stability and nonhyperbolicity. 

We point out that although power-laws such as that given in Eq. ( !22|) are susceptible to 
variations in the initial position {qo,Po), to changes of the trajectory sampling, and to the 
time window in which they are under study, the a deviations from unity are robust structures 
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FIG. 7: Small regular islands found at K = (8.47, 9.35) and their surrounding homoclinic tangles. 
Not only are the islands much smaller than for K = 9.26, but the turnstile fluxes are larger relative 
to the measure of the enclosed area, making the flow much less restricted in the neighborhoods of 
the islands. The measure of the regular region for K = 8.47 is roughly 0.01% of counting all the 
unique iterations of the structure. 

and can be correlated to phase space structures. In Fig. (JHl) we show the variation of a as 
is tuned through a range which has small stable regions. The "normal " exponent a = 1 is 
realized to a very good approximation for most values of K indicating the general absence of 
significantly large, low instability traps. The prominent deviation around K = 9.3 is due to 
the structure of Fig. (jlj) born out of a double saddle node bifurcation around K = 9.2 and 
it persists to near K = 9.45. It grows to its largest size around K = 9.26 where the variance 
decays the most slowly, a ~ 0.16. In terms of the phase space area, the islands occupy 
approximately 0.5% of n at K = 9.26 and 0.2% at K = 9.35 (for which a ^ 0.88). Notice 
robust, though small, deviations from unity at several other values of K. For example, even 
for an island as small as that found at i^' = 8.47 (0.01%, a ~ 0.96), there is a range of K 
values where a is significantly and consistently different from unity; the inset shows a finer 
scan of this region. The fine oscillatory structure that is apparent here could be reflecting 
the fractal behavior of such exponents themselves which in turn are due to bifurcations 
in the map that produce or destroy very tiny islands almost continuously. The Newhouse 
phenomenon wherein homoclinic tangencies produce an infinite number of sinks is replaced 



in conservative scenarios by a large number of elliptic islands 32|. Knowing that there is 



possibly an island somewhere in phase space does not locate it, and a blind search could 
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FIG. 8: The exponent a in the decay of the variance of the FTSE as a function of the kicking 
strength K. The ensemble is the same as that used in Fig. ([6]). A region around K = 8.47 is shown 
at a resolution of 0.001 in K. 

well be futile. Some direction is important. The islands are found more easily by searching 
around regions of the lowest stability exponent, regions that are highlighted in figures such 
as Fig. ([3]), and correspond to zones of low hyperbolicity. Note one last observation, for the 
three islands constructed, K = (9.26, 9.35, 8.47), the degree of a's deviation from unity and 
the measure of the islands are in corresponding order, smallest to largest. 

C. The higher order cumulants 

Similar more complicated dynamical correlation expressions enter for all the higher cu- 
mulants. They naturally weight more the extremes of the probability density, and thus, 
one might expect them to be even more sensitive to stable islands. To be useful, they also 
need a standard time scaling for their correlation independent components as found for the 
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variance. In fact they do have a simple time dependence and it involves generalizations of 
X2. By defining: 



the correlation-independent component of the higher order cumulants are expressible as 
polynomials in these quantities with a single overall power-law decay in time and no other 
parametric dependencies. Details are presented in Appendix B, but note here an interesting 
consequence of these expressions that imply a special feature found in any of the 2D kicked 
rotors. Substitution gives 



which shows that to the extent that dynamical correlations can be ignored, every cumulant 
is independent of the overall kicking strength. The probability density thus has a single 
functional form for all kicking strengths. Therefore, although the Lyapunov exponent grows 
logarithmically with K, the probability density of FTSE retains the same shape and width 
(it just translates with ui as K increases). Assuming the kicking strength is large enough 
that the system is strongly chaotic, any K dependence of the probability density, other than 
where it is centered, is due to variation of dynamical correlations and oscillatory. 

Even stronger dynamical correlations and more sharply defined fluctuations show up in 
the higher cumulants. In Fig. ([9]), the fourth and second cumulants are compared. The 
fourth cumulant displays much sharper features than the second. On the other hand, the 
features are present in the same regions of K, so the connection to stable islands is the same 
as for the variance. We did not carry out the corresponding time scaling calculations as we 
have no reason to doubt that they reflect the stable structures for the same reasons as for 
the variance. 

D. The fluctuations in l/|Det(Mt - /)| 

Since higher cumulants carry information and have increased sensitivity to stable struc- 
tures, a quantity that contains information from all the cumulants has the potential to be 
roughly optimal for its sensitivity to stable structures. Through relations such as Eq. (ITSl) . 
one sees that l/|Det(Mt — /)| fits into this category. Furthermore, it is a quantity of consid- 
erable interest in semiclassical theory, where the square root occurs as the amplitude along 
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FIG. 9: The fourth and second cumulants as a function of K using the j = approximation. The 
sohd hne is the fourth cumulant and the dashed Une is the second. 

classical paths. We therefore define 

^'-T'"( |Det(M-/)| ) 
In a uniformly hyperbolic system, (3oo is the Lyapunov exponent. However, if the hyper- 
bolicity is not uniform, as is generically the case, then f3t is less than the corresponding 
stability exponent. In the presence of stable motion, it must tend toward vanishing. The 
averaging above can be done over any ensemble of phase space trajectories over time t\ we 
use a uniform density over VL. Shown in Fig. (fTOl) is the variation of this quantity across a 
range of K values. Its extreme sensitivity to small islands and sticky regions of phase space 
is quite apparent. This may be expected as the denominator can come arbitrarily close to 
zero in the neighborhood of points undergoing bifurcations. We show two calculations in the 
figure, one of which is an "exact" evaluation based on finding the trace of Mt at large times, 
using a determinant method used initially for periodic orbits |33|, l3J] , and the other is by 



using the product (j = 0) approximation in Eq.([7]). The latter procedure is a stringent test 
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FIG. 10: The quantity f3t as a function of at t = 100. The averaging involved 40, 000 initial 
conditions spread uniformly across phase space. The solid line is the exact trace. The relatively 
smoothly behaving dashed line is what the result would be for a normal variable, i.e. (3t = I'l — if2/2. 
The dotted line derives from the j = approximation. In principle, the approximation's errors are 
much worse for a quantity such as /3t relative to one of the low order cumulants, Vn- 

for this approximation and it is interesting that it picks out almost all of the "hotspots" . It 
is however much faster and easier to implement than the exact method. The smooth line 
shown is what one may expect from a normally distributed variable. 

On comparison with a and with the variance itself, we find that the ranges of K where f3t 
deviates from that expected of a normal random variable, is also reflected in these quantities. 
Note the consistency in all these measures with respect to structures due to stable, if small, 
islands in the phase space. For example, the structures around K = 8.47 reflect those already 
illustrated in Fig. ([8]), while those around K = 8.66 correspond to similar islands. We again 
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expect that higher dimensional generahzations of the approximation j = studied here will 
provide a method of ascertaining higher dimensional, hard to detect, 'weakly' mixed phase 
spaces. 

V. DISCUSSION 

The consideration of finite-time stability exponents has some advantages vis-a-vis finite- 
time Lyapunov exponents. For one, they are more directly related to quantities that ap- 
pear naturally in quantum mechanics. For another, their generalization to many-degree-of- 
freedom systems tend to the Kolmogorov-Sinai entropy in the t —>■ oo limit. There is a local 
approximation just as there is for FTLE, but the time-dependent transient may disappear 
exponentially fast in many cases for FTSE and not FTLE, depending on averaging tech- 
niques employed. Corrections to the local approximation are also relatively simple and the 
first one is sufficient to calculate improved estimates of Lyapunov exponents. We showed 
that this gives a leading correction to In y of exactly (K'^ — 4)^^ for the standard map. In 
fact, although not explicitly mentioned, we explored a wide variety of techniques all of which 
failed to give this correction properly (most vanished or diverged). The corrections them- 
selves are interesting as well. One would guess that they should lead to a Levy-process and 
numerically this is borne out. At least for low-order cumulants in the large kicking strength 
regime, the local approximation along with the first correction gives the full picture, and for 
many quantities the local approximation is more than sufficient. 

For the standard map, the mean (first cumulant) FTSE converges rapidly to the Lyapunov 
exponent. This is understood to be a consequence of Cq tending exponentially quickly to 
zero with time in Eq. f|T71) . The local approximation gives a mean FTSE that is independent 
of direct dynamical correlation contributions. Indirectly, in the mostly chaotic regime, the 
exclusion of chaotic trajectories from small regular regions (islands) would slightly alter the 
ergodic averaging, but never more than the relative measure of those regions. Consistent 
with this picture, other than small deviations, one sees very little effect of small islands on 
the mean FTSE. 

The fluctuations behave quite the opposite. They can be investigated through the creation 
of a probability density for the FTSE. In the absence of small regular regions, the probability 
density takes on a single form, independent of system parameters (such as the kicking 
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strength), which slowly shifts towards a Gaussian form as t — oo of width 7r^/12. Therefore, 
all deviations from this form must be due to dynamical correlations, and in particular, 
correlations introduced by the existence of small stable islands of motion. Since very strong 
K dependence is seen, the higher cumulants are very sensitive sensors of the presence of such 
non-hyperbolicity in the system. As such, they can be used as 'detectors'. An important 
quantity which shows up in semiclassical theories, Det|Mf — in principle, contains 

information carried in all of the cumulants. In fact, it may be the most sensitive detector 
of dynamical correlations induced by tiny regions of non-hyperbolicity in the dynamics. All 
of the techniques and main results discussed have their generalizations to many degree-of- 



freedom systems. Those studies are currently underway [20 1. 
Appendix A 
The relation is 

odd t — 2j 
[r/ ■^J arrng. factors 

Trmqo,po)] = Y.{-iy n [2 - (26) 

j=0 i k 

where [t/2] is the integer part of t/2 and the odd arrangements of the product of t — 2j 
factors is specified just ahead. Also note that if t is even, the final term (the one for which 
3= t/2) is (-l)*/2 X 2. 

There are two key properties required for the above identity. First, there must exist at 
least one column or row in the set of {M(gj)} both of whose elements are the same constants 
for all members of the set. Second, the determinants of the set of {M(gj)} must all be equal 
to the same value (if this value is not equal to unity, then the terms in the equation above 
must be multiplied by its j^^ power) . 

It turns out there is a fairly simple rule for which products can contribute to a particular 
"j-term". Imagine the t integers {i G [0,i(: — 1]} uniformly spaced on a ring. Those are the 
"locations" that each factor in a product can occupy, so-to-speak, determined by its own 
subscript. For a given time t with t — 2j factors in the product, each arrangement that has 
exclusively odd-unit-spacings between occupied locations contributes once to that j-term. 
All cyclically connected arrangements are automatically included in this construction, and 
thus the cyclic invariance of the trace is easy to see. The number of such odd spacing 
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arrangements Ma turns out to be 



Ma 



t ft- 3 
t-3 V J 



(27) 



which is related to the summation and recursion relation 

1 d-'' 



[i/2] 
m=j 



(1 + v^) +(l-v^)' 



t 

2m 



t-3 



x=l 

t-3 
3 



(28) 



The combinatorial are related to an expansion of the functional form of eigenvalues raised 
to various powers of a two-by-two matrix. The j = term is simply the product of all 
factors 

Tr[M,(go,Po)],=o = H " ^"(^O] (29) 

One could also label this term {1*} indicating that there are t unit spacings. Similarly, the 
3 = 1 term is always the cyclic sum 3} of which there are t terms. Thus, 

t-i 

Tr[Mt(go,Po)]j=o+i 



\{[2-V"{q,)]x 



cyclic 



1 



i=0 



[2 - V"(g,)] [2 - V"{q,^,] 



(30) 



where the last term has indices {t — 1,0). Pulling the product of all factors out front 
emphasizes that if the {V"{qi)} are generally much larger than 2, then for times not too long 
the j > terms can be thought of as corrections, each successive one roughly proportional 
to t < (y"Y weaker that the previous term. In some sense, the first non-trivial j- 
term comes for t = 6 and j = 2. One expects 9 double products of factors. The allowed 
odd-spacing arrangements are {1,5} and {3,3}. There are 6 terms related to the cyclic 
permutations of the i{i + 1) product ({1,5}) and only 3 terms for the i{i + 3) product 
({3, 3}) since it repeats itself after half a cycle, and thus 9, the correct count of terms, is 
consistent with Eq. ( 1271) . The above expression is the expansion of the determinant given 
by Bountis and Helleman 33|| nearly 30 years ago. Greene developed its small-i^ expansion 
in his study of the behavior of KAM (Kolmogorov-Arnol'd-Moser) surfaces 34 1. 



Appendix B 



It is more convenient to focus on the cumulants than the central moments as it simplifies 
the time dependence. There is also the advantage of interpretation of the probability density. 
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FIG. 11: The first four cumulants as a function of time. The cumulants are multiplied by the 
time factors needed to make them independent of time according to the expressions leaving out 
dynamical correlations. The kicking strength is 8.9633781407, which gives the prediction that the 
Lyapunov exponent (first cumulant) is equal to 1.513. 

For a Gaussian density, only the first two cumulants are non-zero. For other densities, the 
higher cumulants describe shape deviations. There is a third advantage, which comes into 
play when considering the relationship with quantum mechanics, but we do not elaborate 
here. If the dynamical correlations vanish or are small enough to be ignored, then the 
expression for the cumulants can be simplified. With these quantities, the cumulants have 
the same polynomial form, P„, as with the central moments and a very simple dependence 
on t. The result is 



Unit) = t^-^P, 



(31) 
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For example, 

Mt) = ^ (^^4 - 3A^|) (32) 

An amusing, known consequence of the time dependence in Eq. fl31l) is that there is a 
sense in which the distribution collapses to a (5-function density as t — oo as expected 
(every trajectory has the same Lyapunov exponent), and the time-dependence of the reduced 
cumulants goes as 

-nit) = ^ oc t-"/^ (33) 

^2 it) 



which for ?7,>3^0ast^oo. In other words, the density slowly approaches a Gaussian 
form. Nevertheless, the 5-function limit is approached slowly enough that the width of 
the probability density for |Tr(Mt)| grows without bound as t increases, and it is not well 
approximated by a lognormal density. An example of the time dependence of the first four 
cumulants for a kicked rotor can be seen in Fig. ffTTj) . It illustrates that the time dependence 
is to within fluctuations as predicted for K values at which dynamical correlations are not 
observed. 
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